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g ■ ABSTRACT 
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ZJ \ The propagation and evolution of planet-generated density waves in protoplanetary 

' disks is considered. The evolution of waves, leading to the shock formation and wake 

P\J ' dissipation, is followed in the weakly nonlinear regime. The local approach of Goodman 

&: Rafikov (2001) is extended to include the effects of surface density and temperature 
variations in the disk as well as the disk cylindrical geometry and nonuniform shear. 
, Wave damping due to shocks is demonstrated to be a nonlocal process spanning a 

I significant fraction of the disk. Torques induced by the planet could be significant drivers 

I of disk evolution on timescales ~ 10^ — 10^ yr even in the absence of strong background 

viscosity. A global prescription for angular momentum deposition is developed which 
could be incorporated into the study of gap formation in a gaseous disk around the 
. planet. 

Q . Subject headings: planets and satellites: general — solar system: formation — (stars:) 

planetary systems 



(N 



1. Introduction. 

Tidal disk-companion interactions are important in a variety of astrophysical contexts ranging 
from the formation and evolution of protoplanetary systems to the origin of galactic spiral struc- 
ture. Gravitational interaction between the disk and companion generates density waves in the 
disk (gaseous, stellar or particulate) which carry angular momentum. This angular momentum is 
eventually deposited elsewhere in the disk, leading to the evolution of the disk as a whole. 

In the case of protoplanetary systems disk-planet interactions not only cause the migration of 
the planet (Goldreich & Tremaine 1980, hereafter GT80; Ward 1997) but can also drive noticeable 
evolution of the disk itself (Larson 1989; Goodman & Rafikov 2001, hereafter GROl). However, to 
change the state of the disk and orbit of the planet in these systems it is necessary to somehow 
transfer angular momentum from the density waves launched by the perturber to the disk material 
and this can only be accomplished by virtue of some damping process (Goldreich & Nicholson 
1989). 
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Various mechanisms have been envisaged as possible sources for such damping. The most 
popular linear ones are viscosity in the disk (Takeuchi et al. 1996) and radiative damping of the 
tidal perturbations (Cassen & Woolum 1996). Viscosity can dissipate tidal perturbations on scales 
smaller than the typical disk sizes only if it is large enough, a > 10~^ (Takeuchi et al. 1996; GROl), 
but it is hard to identify a strong source of viscosity in protoplanetary disks. The most probable 
viscous mechanism in hot accretion disks - magnetohydrodynamic (MHD) turbulence driven by the 
magnetorotational instability (Velikhov 1959; Balbus Sz Hawley 1998) - probably does not operate 
in protoplanetary disk environments: the gas is too cold and weakly ionized throughout most of the 
disk (Jin 1996; Hawley & Stone 1998). Convection was put forward as another possible source of 
viscosity (Lin & Papaloizou 1980) but analytical arguments and numerical simulations cast serious 
doubt on the ability of this mechanism to produce outward angular momentum transport (Ryu & 
Goodman 1992; Balbus et al. 1996). The efficiency of radiative damping in protoplanetary disks 
is strongly reduced by dust opacity (Henning & Stognienko 1996), which leads to very high optical 
depths (r > 10^) and implies very small radiative losses from the disk surface. Waves could also be 
damped by radiative transfer in the plane of the disk but this turns out not to be very important 
either (GROl). 

Nonlinear dissipation, namely shock formation and its consequent damping, seems to represent 
a more efficient and almost inevitable process for transferring the angular momentum of the wave 
to the disk fluid. There are two reasons for this. First, the differential rotation of the background 
fluid causes the wavelength of the tidal perturbation to decrease as it travels away from the planet 
(Goldrcich & Tremaine 1979). Second, the amplitude of the planet-induced wake is growing, at 
least in the planetary vicinity, as a consequence of the conservation of the angular momentum 
flux. These processes working together can make the wave shock very rapidly if the initial wave 
amplitude is significant (i.e. if the perturber is massive enough). 

Goodman & Rafikov (GROl) have performed a detailed study of nonlinear evolution of planet- 
induced density waves in two-dimensional disks. Using the shearing sheet approximation they have 
demonstrated that for small enough planetary mass it is possible to separate the wake evolution into 
two distinct stages: linear generation, which takes about {1 — 2) x h from the planet to complete 
(here h = c/Q is a typical scale length, which cqTials disk thickness in three dimensions, c is a 
sound speed and 17 is a disk angular frequency), and then nonlinear evolution of the wake causing 
it to shock after travelling several h from the perturber (of course, depending on the mass of the 
planet). After the shock is formed, it damps transferring its angular momentum to the mean fiow 
and leading to disk evolution. In favorable conditions the damping of the density waves can make 
the disk evolve on timescales comparable with those derived from observations [10^ — 10^ yr, see 
Hartmann et al. (1998)]. 

The considerations of GROl were in a certain sense local because of the shearing sheet ap- 
proximation. It was shown that this approach is good for the description of the shock formation 

but probably not so accurate for studying the shock dissipation: in the shearing sheet geometry 
shock damping proceeds slowly and at some point the background fluid velocity can no longer be 



-3- 



represented by just a uniform shear. The analysis of GROl also assumed a disk with constant 
background surface density and sound speed and did not take into account the effects of the disk's 
polar geometry on the evolution of the density wave amplitude. These approximations naturally 
lead to a picture in which the wake itself and its damping pattern are symmetric on both sides of 
the planet. 

The main purpose of this paper is to extend the analysis of GROl by including the effects of 
radial surface density and sound speed variations in the disk on the behavior of the weak tidal 
disturbances generated by a low-mass planet incapable of opening a gap in the disk. We consider a 
Keplerian rotation law (not a linearly sheared background flow) and a polar geometry to include self- 
consistently important ingredients needed to provide a global picture of the nonlinear evolution of 
the density waves in non self-gravitating disks. We restrict our attention to purely two-dimensional 
disks, thus completely disregarding vertical motions and related phenomena, such as wave-action 
channeling in the vertical direction (Lin et al. 1990; Lubow & Ogilvie 1998). We believe that this 
is a good approximation in passive, externally irradiated protoplanetary disks with high optical 
depths which should have almost isothermal vertical structure (Chiang k, Goldreich 1997). 

The paper is structured as follows. In §2 we study the full system of fluid equations in cylin- 
drical geometry. We provide in §2.1 simple scaling laws for the behavior of the wake amplitude 
and shock formation distance in the quazilinear approximation from rather general qualitative con- 
siderations. We then confirm these estimates by an accurate quantitative analysis in §2.2. Wake 
properties are studied in §3, in particular for the case of disks with power-law surface density and 
sound speed radial dependencies (§3.2). Applications of our results are discussed in §4. 



We consider a system consisting of a gaseous non self-gravitating disk rotating in the unper- 
turbed potential C/*(r) (r is the distance from the central object) and a planet of mass Mp located 
at a distance Vp from the center. The disk is assumed to be geometrically thin and its unperturbed 

surface density Eq and sound speed cq are both taken to be functions of r. Disk scale height in three 
dimensions is h = cq/VL. We will always denote by the subscript "p" various quantities evaluated 
at the position of the planet. 

The equations of motion and continuity for the two-dimensional disk read as always 



Here v and S are the fluid velocity and surface density, P is the pressure, and the potential 



2. Density wave structure. 
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consists of three contributions: the potentials of the central star and perturbing planet (we assume 
a Keplerian potential = — GMc/r, where Mc is a mass of the central star), and the indirect 
potential Ui, which can always be neglected here (because disk mass is much smaller than Mc). 

We assume that the pressure P is related to the instantaneous value of the surface density E 
by the locally polytropic law with a specific index 7: 



Poir) 



1 1 



So(r)J 
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Then the perturbed sound velocity is given by the usual expression 
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meaning that Po{r) = Eo(r)co(r)/7. Clearly, Po{r) and co(r) are unperturbed values of pressure 
and sound speed. The equation of state given by (4) does not force co(r) to be related to Eo(r) and 
this gives us additional flexibility in applications. The entropy of the disk fluid now varies with r, 
contrary to the equation of state with a radially-independent polytropic constant. 



2.1. General considerations. 

Using basic physical principles such as conservation of angular momentum flux it is possible 

to analyze the behavior of the wake in the quazilinear regime and determine when the wave shocks 
on a qualitative level. The simple scaling laws obtained in this way will later be confirmed by more 
rigorous analysis of the complete system of the fluid equations. 

Let us first concentrate on linear wake propagation. One can easily show that the solution 
of equations (l)-(2) in the linearized form with the equation of state (4) yields, using the WKB 
approximation, 

{n - Qpf = K^ + k^c^, (6) 



and 

d_ 
dr 



3 

-(S-Eo)^ 



for non self-gravitating disks. Here m is the azimuthal wavenumbcr, k is the radial wavenumber of 
the perturbation, VLp is the perturbation pattern angular frequency, (E — Eo)m is the m-th harmonic 
of the surface density perturbation, and 

0^ = lf^ + -L^, n\r)=AB{rmr), (8) 
r dr En dr 



where 



(9) 
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Equation (6) is the usual dispersion relation for small perturbations in the thin differentially- 
rotating disk and it shows that density waves behave basically like sound waves after propagating 
several scale lengths h from the perturber. Both equations (6) and (7) coincide with their analogs 
for the disks with constant entropy (Goldreich & Tremaine 1979). One can also demonstrate that 
the total angular momentum flux across a cylinder of radius r carried by the m-th harmonic of 
surface density perturbation is given by 

[which is the same as the analogous expression of Goldreich & Tremaine (1979)], which, combined 
with equation(7), implies that fj is conserved. 

If the nonlinear effects are fully neglected and fj is strictly conserved as demonstrated above, 
then it follows from (10) that the magnitude of the surface density perturbations scales with other 
flow variables as 

(S - So)^ oc (11) 

This equation demonstrates how the amplitude of the wave varies in the linear regime and we will 
return to it in §2.2. 

In reality, even if the perturbation is small, different points of the wake profile have different 
propagation velocities so that the waveform constantly distorts. This nonlinear evolution leads to 
shock formation. In a shock the angular momentum of the density wave gets transferred to the 
mean flow so that scaling provided by (11) breaks down after the shock is formed. We can estimate 
when this happens. 

Let us consider a part of the profile which initially has a phase = (po. It evolves according to 



dip 
dr 



^0 A(r)c(r) 



where Sc{r) is the perturbation of the propagation velocity c, which is the sound velocity in our 
case (it is different for different ipo and this is responsible for the profile distortion), and A(r) is the 
current wavelength. 

Integrating (12) with respect to r one gets that 

= (^0 + 27r / -TpYY^dr (13) 
J X(r)c{r) 

ro 

(ro corresponds to the point where the wave is launched). 

The wave shocks when its profile acquires an infinite slope after travelling some distance rgh- 
d(p/dipo = at r = Vgh for some specific tpo. This gives us [using (13)] the condition for Vgh in the 



following form: 



rsh 



I 



ro 



k{r)5c{r) 
c(r) 



dr = const. 



(14) 



Prom the dispersion relation (6) we get that A; oc ($7 — ^p)/c; since Sc/c oc (5S/S, we can use the 
scaling provided by equation (11) to obtain finally the shocking condition in the form 



rsh 

I 

ro 



0^0 



dr = Cfj^''' oc M;\ 



(15) 



where C is some constant. 



Scaling laws (11) and (15) immediately provide important information about the wake propa- 
gation in the linear regime and conditions needed for the shock formation. 



2.2. Basic equations. 

In this section we study the behavior of the wake by solving the fluid equations for small 
perturbations in a weakly nonlinear regime. We work in the polar coordinate system rotating with 

the angular velocity of the pcrturbcr Jlp = ^J {1 / rp){dUi, / dr)\rj,. In this coordinate frame the flow 
is stationary (time-independent). As always, we take the r-axis to be directed out from the central 
body and the ^-axis to be directed in a prograde sense. 

In this coordinate system the equations of motion in the r and ^ directions are (Landau h 
Lifshitz 1959) 

VrdrVr + —d^Vr ^ = — — drP + 20pf<A -|- O^r — drU^ (16) 

VrdrV^ + '^d^V^ + = -^9^P - "^^pVr - ^d^U, (17) 

and the continuity equation is 

-dr i^rvr) + -d^ i^v^) = 0. (18) 
r r 

Here Vr and v^f, are the fluid velocities in the r and (f) directions respectively. Equations (16) and 
(17) include Coriolis and centrifugal forces. 

GROl have demonstrated that the system (16)-(18) simplifies significantly if the mass of the 
perturber Mp is smaller than a characteristic mass M\ given by 
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where A{r) = {r/2)dQ/dr is the Oort's A constant which is related to B{r) hj B = Q + A 

[equation(9)]. Then the disk could be split into two distinct regions: an excitation region, within 
several scale lengths hp = h{rp) = CQ{rp)/rt{rp) from the planet, where one can neglect nonlinear 
effects, and a propagation region beyond several hp from the planet, where the planetary potential 
is negligible and one can study wake evolution caused by nonlinear effects. 

In the first region numerous Lindblad resonances tidally excite density waves corresponding 
to different azimuthal harmonics and provide individual Fourier contributions to the angular mo- 
mentum flux (GT80). Most of the angular momentum comes from the resonances with azimuthal 
wavenumbers m ~ rp/hp S> 1 which lie close to the planet — at distances ~ hp from it. Harmonics 
with smaller m are weaker because the tidal excitation they experience is reduced by their larger 
distance from the planet. Fourier components with m higher than Vp/hp are strongly damped 
because of the so-called "torque cutoff" (GT80). Its nature can be qualitatively understood as 
follows: near the planet, less than {2/3)hp from it (in a Keplerian disk), the background fluid flow 
is subsonic which does not allow a stationary perturber to excite sound waves (Landau & Lifshitz 
1959). This cutoff strongly reduces the amount of the angular momentum flux carried by the 
corresponding harmonics. 

GROl studied the linear wake formation process in the shearing sheet approximation taking 
into account not only the contributions of individual resonaces to the torque, but also the phases 
of the Fourier harmonics of the surface density perturbation which allowed them to obtain the 
shape of the wake in the linear theory. Here we will simply use their results in our more general 
case because wake generation is essentially a local process, spanning only a few hp from the planet, 
where the shearing sheet approximation gives a good representation of the disk velocity profile and 
So and Cq may be assumed constant. The solution of this problem in our more global setting with 
varying Eq and cq will only be different from the previously obtained local solution by factors of 
the order of hp/vp <^ 1 which is of no interest for us here. 

Thus, we can proceed immediately to study the wake propagation region. To do this we use 
a simple extension of conventional perturbation theory capable of including a weak nonlinearity of 
the wave. We assume that Vr and Vfj, are given by 

Vr = u, V(f, = vo{r) + v with vo{r) = r [Q{r) — ilp] , (20) 

where u and v are the velocity perturbations and we take u,v <C vq since the shock is assumed 
to be weak and we are always several scale lengths away from the planet. We will often write for 
brevity AJ7 = D,[r) — Qp. 

Substituting (20) into (16)-(18) one obtains the following perturbation equations: 

/I 1 \ v"^ 

And^u + udru -2^v = -[ -drP - ^drPo j + — > (21) 

1 vu 

Andj^v + udrv + 2Bu = — -d^P , (22) 

r 
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5j 5j 

Ands^ + udr'S + = -9<aS dsv. (23) 

We have everywhere neglected v in comparison with Af^r and made use of the fact that c^^Pq = 0. 
In equations (21)- (23) terms quadratic in u and v are subdominant; however we wih keep ud^u 
terms because they are the strongest drivers of nonhnear evolution. Also, we assume that v <^u 
and dtf, <^ rdr as a consequence of tight-winding approximation. These assumptions are checked in 
Appendix. 

By introducing a new radial coordinate given by 

r 

i = J [n{r) - Qp] dr (24) 

rp 

this system is transformed into 

11 2Q v"^ 

d,u + ud,u + -d,P - -d,Po = ^v + ^^, (25) 

+ uO,Y + = - - ^a,.. (27) 

The l.h.s. of (25) and (27) are similar to the usual system of equations describing isentropic 
one-dimensional gas motion (Landau & Lifshitz 1959), which possesses two invariants conserved 
on characteristics — Riemann invariants R±. In our case the nonzero r.h.s. of (25) and (27) cause 
these invariants not to be conserved exactly, but one can still use them in a slightly modified way. 
We extract from (25) and (27) two equations of evolution of the Riemann invariants: 



[d^ + {u±c)d^]R± 



-d,P - -d,P, - cd,^^ ± cu— T ud,^^ 
20, cu cv ^ , ^ c „ 



where 



R^ = u±^^, (29) 
7-1 

and it is always assumed that c = c(r) = c{S,). 

As noted in GROl, the characteristic C+ crosses the wake profile while the characteristic 
C_ follows it, that is C_ is always in the perturbed region while C_)_ mostly passes through the 
unperturbed fluid and only for short periods of time enters the perturbed region. One can assume 
the changes in due to these crossings to be small and take R+ to be constant (we will later 
comment on the effect of the multiple crossings) and equal to its value in the unperturbed region. 
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There, i?+ = 2co/(7 ~ 1) everywhere. Then, from the conservation of one obtains that u = 
2(co - c)/(7 - 1) and = 2(co - 2c)/(7 - 1). 

Using these results it is demonstrated in Appendix that the equation of evolution of i?_ can 
be reduced to an inviscid Burger's equation 



where 



7 + lS-So , , 

X = -i; — gin, 



t = 



Iv J co{r')g{r 



ip 



(j) + sign(r — r. 



r 

'^1 



and g{r) = 



0(rO - Op 
co(r') 



dr' 



2i/4 



(30) 

(31) 
(32) 

(33) 
(34) 



Here Ip = Cp/\2A{rp)\ = {2/3)hp is a Mach-1 length (distance from the planet where the Keplerian 
shear makes fluid velocity equal to Cp). 



In the immediate vicinity of the planet (but still several hp from it) definitions (31)-(34) reduce 



to 



X 



7 + 1 



23/4 



-1/2 



23/4 



5/2 



y X ■ / \ 



(35) 



where x = r — Vp and y = rp(f). These equations coincide with analogous expressions in GROl. 

The tidal perturbation launched by the planet follows a nearly parabolic path in the immediate 
vicinity of the planet where the shear can be assumed constant. Further from the perturber where 
the shear is no longer uniform, the density wave has a spiral shape. Its pattern is described by the 
equation 



bo — sign(r — Vp 



'I 



O(r') - ^p 
co(r') 



dr' 



(36) 



Depending on the conditions in the disk this spiral can wind up several times around the center 
before the wave damps. One can see from (36) that the perturbation is indeed tightly wound if 
|r — Tpl ^ hp (this could, however, be violated further away from the planet for some profiles of Sq 
and Co, see §3.2). 
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Now, if we disregard the nonlinear evolution entirely, it follows from equation (30) that x 
fkn) [/(^) is an arbitrary function of the argument x], or 



(S - So) 



lin 



cog{r) 



firi) 



rcr, 



1/2 



/ <^ + sii 



ign(r-rp) j 



co(r') 



(37) 



This means that in the linear regime, the waveform propagates along the wake (whose location is 
given by the condition rj = const) and its amplitude scales with the distance from the planet in 
complete agreement with equation (11). 



3. Nonlinear evolution of the wake. 

Burger's equation (30) is probably the simplest partial differential equation able to exhibit 
a shock formation phenomenon. For the reasons outlined in §2.1 the profile of the density wave 
produced by the linear generation mechanism is constantly distorted in the course of its propagation 
away from the planet, so that finally the waveform breaks to become double- valued. This implies 
that a shock must be formed at this point. The distance which the density wave travels before it 
shocks depends on the initial shape and the amplitude of the wake. 

To study shock formation and propagation quantitatively we need to solve equation (30) with 
the initial condition given by the solution of the linear wake generation problem. As we have 
mentioned in §2.2, since all our variables reduce to the corresponding variables of GROl in the limit 
r ^ Tp, the whole linear generation problem reduces to the one solved before — wake excitation in 
the shearing sheet approximation. Thus, we might use the solution for the wake shape calculated 
in GROl as the initial condition in our more general case: 

x{Mp,t = to,r,) = ^Jo{v) (38) 

[Ml is defined in equation (19)]. Here the initial profile x(Mp,t = to,rf) is taken not at t = but 

at some to > because in the linear regime wake has to propagate some distance from the planet 
before it fully forms. Thus is a matching boundary where one switches from linear to weakly 
nonlinear regime, and following GROl we take to ~ 1-89 corresponding to the distance from the 
planet \xo\ = 2lp = (4/3) Zip. One can usually neglect to once the wave has travelled several hp away 
from the planet. 

The function fo{ri) represents the shape of the initial profile at the moment t = to for Mp = 
Ml [see equation (19)]. It was calculated by GROl in the shearing sheet approximation and is 
reproduced for convenience in Fig. la. The factor Mp/Mi rescales the amplitude of the wake for 
an arbitrary mass of the planet. 

In terms of the variables x, r), t our nonlinear problem (30) is identical to the one studied in 
the shearing sheet approximation [including the initial condition (38)]. This means that the result 
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of GROl for x{Mi,t,r]) in terms of variables x,r7,t is also a solution of our more general problem 

in terms of these variables. One needs only to express them in terms of r and <p using equations 
(31)-(34) and rescale to an arbitrary mass Mp in the following way: suppose that xii^ — toiV) = 
x(Mi,i — to,r]) is a solution of the equation (30) with the boundary condition (38) and Mp = Mi 
(calculated in GROl). Then one can easily see that the solution for an arbitrary Mp can be written 
as 

xiMp, t - to, v) = ^Xi (^(i - to), ■ (39) 

This reduction to the previously studied case is a very convenient feature of the analysis 
because it allows one to use all results previously obtained by GROl in our more general setting. 
For instance, it was found before that the shock develops at 

Ml 

tsh = to + 0-79-^ (40) 

for the initial profile given by the function /o(??) (see Fig. la). Based on what we have previously 
said, the shock has to form at the same value of tgh in our global case, and the corresponding radial 
distance from the planet rgh — fp at which this happens can be found using equations (32) Sz (34). 
For a fixed Mp this procedure leads to the condition (15) found before on the basis of qualitative 
considerations. 

After the wave shocks, it starts transferring its angular momentum to the disk mean flow which 
leads to the damping of the wave amplitude. The initial profile shown in Fig. la has both positive 
and negative parts, thus it is destined to evolve asymptotically into an "N-wavc" profile (Landau 
& Lifshitz 1959; Whitham 1974). In this regime the amplitude of x decays as t^^^^, whereas it was 
constant near the planet before forming a shock. The width of the spiral perturbation, which was 
^ hp initially, grows as t^/^ in the "N-wave" stage. 



3.1. Angular momentum flux. 

Summing up all the angular momentum flux contributions given by equation (10) and applying 
Parseval's theorem one can show that the total angular momentum flux of the density wave fj is 
given by: 

TT 
—IT 

Using definitions (31) and (34) we obtain that in terms of variables t and rj 

23/2^3^ 

*(Mp,t) = |x'(Mp,t,>,)<i>,= (^) *(m,,-^(), (42) 
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1 




Fig. 1. — (top) Shape of the wake profile at the moment t = to ^ 1-89 calculated in the linear 
theory. This profile is used as initial condition for the calculation of the nonlinear evolution of the 
wake, (bottom) Behavior of the dimensionless angular momentum flux $(t) carried by the waves in 
the presence of shock damping. The calculation is done for Mp = Mi . The dashed line has a slope 
of t~^/^ and is drawn to illustrate the asymptotic behavior of ^{t). 
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exactly like in GROl (we neglected to here). Prom equation (30) one can easily see that without 
shocks $(Mp,t) = const and angular momentum flux is conserved. 

After the shock forms $ starts to decrease. Asymptotically, in the "N-wave" regime it falls like 
$(Mp,t) oc The behavior of $ as a function of t for M.p = Mi was calculated in GROl and is 

shown in Fig. lb. What happens with the angular momentum flux of the wave in the disk depends 
on how t is related to the distance travelled by the wave. If one assumes a uniform shear, then 
according to equation (35) t oc |r — Tpl^/^ and at inflnity the wave damps completely. Real disks 
are different — they are always finite, dissipation can be not complete, and the damping pattern is 
asymmetric between the inner and outer disks. We will demonstrate later that for some background 
surface density and sound speed profiles t{r = oo) or t{r = 0) are finite. This means that a density 
wave that has managed to propagate to the outermost parts of the disk, or to its center, still 
carries some undamped angular momentum flux. It could even happen that t{r = oo) < tgh or 
t(r = 0) < tsh, meaning that the wave does not shock before it reaches the outer or inner edge of 
the disk. The wave could then be reflected and part of its angular momentum could be carried 
back to the planet. This might have important consequences for planetary migration and we will 
dwell upon this more in §4. 



3.2. Power-law disks 



We now consider the case of power-law disks, i.e. Sq and cq are assumed to be some powers 
(usually negative) of r, to illustrate the results obtained before and test some of our assumptions. 
We consider a Keplerian disk with U^, oc and take Sq = Sp(r/rp)~'^ and cq = Cp{r /rp)~^ , 
(5 > 0,z^ > 0. 



Using definitions (31), (32), and (34) we find that 

r/rp 



t = 



and the wake equation 



5/2 



25/4 yh^j 



g3/2 _ ^Z/2^{bv+5)/2-ll/i^^ 



(43) 



.Tp \{r/rpr-^l'' {r/rpr+^ 
(j) = (f)o-sign{r -rp)-^ ^ ^ _ v / z'/ 

fin 



z/- 1/2 



1^+1 (2i/-l)(i/ + l) 



(44) 



Using equation (43) let us consider first the inner (r < r^) part of the disk. One can easily see 
that t^ooasr— >0 if p = 5z/ + (5< Per = 7/2. This means that the density waves launched by the 
planet at Vp in the inner disk span the whole range of t and thus damp completely upon reaching 
the center. Of course, real disks always have an inner cavity, presumably formed by the protostcllar 
magnetospheric activity, where the gas is absent, so we only have to consider wave propagation up 
to this inner edge. Depending on the location of the planet, this inner disk boundary may be very 
close to the center so that we will consider wave propagation down to r = 0. 
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r/rp r/ 

Fig. 2. — Dependence of the dimensionless variable t [multiplied by {hp/vpf'/'^] upon the distance 
from the central body r for several values of the power-law index p = bv + 6 (see the beginning of 
§3.2). Left panel is for r < r^, right one is for r > Vp. Curves are labelled by the corresponding 
values of p. Notice a universal behavior of t{r) near r = Vp, where the shearing sheet approximation 
is valid. In the inner part of the disk t{r) is increasing to infinity as r ^ (which implies a complete 
dissipation of the shock) only for p < pcr = 7/2. In the outer disk wave shocks and damps for all 
values of p as r ^ oo. 
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For p> Per however, r = corresponds to a finite t, and if t(0) < tgh the wave does not shock 
at all as r decreases to 0. Indeed, for p = 5i/ + 5 > pcr- = 7/2 one obtains that 



m 



1 



2V4 \h^) 



5/2 



B 



5z^ + (5 



7 5 
6' 2 
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where B(a,P) is a beta-function. Using equation (40) and neglecting to in it we get that if 



Mp<M^ = MiX 0.94 ( ^ 



5/2 



5u + S 



7 5 
6' 2 



(46) 



then the wave does not shock on the first passage to the disk center. For this to happen Mp should 
be sufficiently small: assuming that B{a,(3) ~ 1, Vp/hp > 1 one needs Mp < Miikp/rpf/"^ for 
this to occur, li p > p^j. but condition (46) is not fulfilled (if the planet is massive enough), the 
wave shocks as it travels towards the disk center, but it does not damp completely upon reaching 
it. Part of the angular momentum flux could be transferred back to the planet in this case. 



One can easily see that in the outer disk this problem never arises: for 5 > Q,u > Q the wave 
always shocks as it propagates outwards, if we forget about the outer boundary of the disk. In 
Table 1 we summarize final outcomes of wave propagation in the inner and outer parts of the disk 
for different values of p. In Fig. 2 we plot the behavior of t{r) for different values of p in the 
inner and outer parts of the disk using equation (43). In the minimum mass Solar nebula (MMSN) 
5 = 3/2 and u = 1/4 (Hayashi 1981) meaning that p = 11/4 < pcr- Thus, in MMSN tidal waves 
always shock and are damped on both sides of the disk. 

It is interesting to see how the shock damping is distributed in the disk. Using the asymptotic 
behavior of $(i) for large t we find that 

/.«(^)'^, (47) 

as f ^ 0. As p varies from to pcr the power law index here changes from 7/8 to (it is equal 
to 3/16 for MMSN), which implies that for some Sq and cq profiles quite a lot of the residual 



Table 1. Wave behavior in the infinite disk for different values of p > 0. 



Outcome of the wave propagation 

Wave shocks, damps completely 

Wave shocks, but damps not completely 

Wave does not shock and does not damp 



Inner disk Outer disk 



< p < Pen any Mp any p > 0, any Mp 

p > Per, Mp>M^ — 
p > Per, Mp<M^ — 



angular momentum is transferred to the disk close to its center. This might significantly enhance 
the accretion rate there. In the outer disk 



for r ^ cx). For shallow profiles of Sq and cq the nonlocality of the damping could be important 
in the outer disk too. The behavior of the dimensionless angular momentum flux $(r) is shown in 
Fig. 3 for different values of p. 

Let us now test the validity of some of our assumptions which were used in the analysis 
presented in §2.2. The tight-winding approximation is one of them. From the equation (36) we 
find that 



where 9 is an angle between the radial direction and the tangent of the spiral, and the last equality in 
(49) is for a Keplerian power-law disk. In the immediate vicinity of the planet, for |r — rp| ~ hp, this 
angle is small and spiral pattern is not tightly wound. But this is the region of the wake generation 
where the free nonlinear propagation approach is not applicable anyway. After travelling several 
hp from the planet, the wake becomes tightly wrapped by the Keplerian shear if the disk thickness 
is small {vp/hp ^ 1). In the outer disk asymptotically ^^7r/2asr— >oo (spiral winds up). 
This happens because in the frame rotating with the planet outer parts of the disk move with 
the angular speed —Up giving rise to a large linear velocity, while cq decreases with growing r. 
Inside the planet's location tan^ pa {rp/hp){r /rpY~^/'^ as r ^ 0. Thus, for v > 1/2 the spiral 
pattern unwinds in the inner regions of the disk and tight-winding approximation might become 
inapplicable. However, if Vp/hp is large enough, wave could still reach the inner disk edge before 
this effect becomes important (in fact, spiral pattern unwinds only if /i ~ r which is usually not 
the case in the inner part of the disk). 

Geometrical effects may also be of some importance. As we mentioned before, in the asymptotic 
regime the wake width increases as t^/'^. At the same time, if the pattern of the wake winds up, the 
distance between the consecutive wavecrests (at a fixed polar angle 0) decreases. At some point the 
"rear" shock front of the "N-wave" profile comes so close to the "forward" shock front of the profile 
lagging by 27r in (j) that our approximation of almost constant Riemann invariant becomes poor 
because the change of during the shock crossing gets comparable to the change of R- following 
the shock (see the discussion in §2.2). 

From Fig. 2 of GROl one can find that the width of the shock in the rj coordinate in the 
asymptotic region t — > oo is Ary ~ 2.3t^/^. Obviously the "front touching" phenomenon occurs 
when {lp/rp)Ar] > 2tt [see equations (AlO) and (A15)], that is when t(hp/rpY'/'^ > V7{hp/rpY/'^ . 
This particular form of the condition is used because it allows direct comparison with Fig. 2. One 
can see from this Figure that for rp/hp = 20 and any p > one can propagate as far as r/r^ ~ 0.3 
in the inner disk and r/r^ > 4 in the outer, and still not encounter this "front touching". In colder 




(48) 




(49) 
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Fig. 3. — Dependence of the dimensionless angular momentum flux of the planet-generated density 
wave $(r) upon the distance from the central body r for several values of the power-law index 
p = 5iy + 6. Left panel is for the inner part of the disk, right is for the outer one (notice that on 
the left panel horizontal scale is linear). Calculation assumes that Mp = 0.2 Mi [see equation (19)] 
and ratio rp/hp = 20. Different curves are labelled by the corresponding values of p. Notice that 
in the case p = 4 > pcr angular momentum flux is nonzero at the disk center. 
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disks with higher r^/hp this hmitation becomes more stringent but even in this case there is a 
significant region of applicabihty of our analysis not too far from the planet where i?+ could be 
assumed approximately constant. 

We also need to mention that our treatment of the nonlinear evolution essentially neglected 
linear dispersion of the wave profile. Near the planet wave dispersion is strongest but it is properly 
taken there into account in the linear calculations of GROl. Further away from the planet, according 
to equation (6), dispersion rapidly becomes less important in comparison with nonlinear steepening 
(if Mp is not very small) and can be ignored. 

4. Discussion and applications. 

It is quite possible that in the presence of vertical temperature gradients the channeling of 
the wave action into the vertical direction can damp density waves in the disk atmosphere more 
efficiently than we have found here (Lin et al. 1990; Lubow k, Ogilvie 1998). It seems reasonable 
however that in passive disks heated by their central stars thermal stratification in the z-direction 
must be small because of their high optical depth. This strongly diminishes the wave action 
channeling into the disk atmosphere (Ogilvie &; Lubow 1999) and leads to almost two-dimensional 
picture of the wave propagation in the disk supporting the validity of our consideration. 

Throughout our analysis we assumed the wave nonlinearity to be weak, meaning that (S — 
So) /So is small. Prom equation (38) we see that if Mp > Mi the wake is nonlinear from the very 
beginning and it shocks immediately [see equation (40)] . This means that the separation of the disk 
into two distinct regions where one can neglect either planetary torques or nonlinearity of the wave 
does not exist. Also for Mp > Mi a gap in the disk can form around the planet (Lin & Papaloizou 
1993). Thus, as we have mentioned before in §2.2, our analysis is applicable only for small-mass 
planets: Mp < Mi. 

We have seen in §3.2 that wave damping can be a nonlocal process and that part of the 
wave action could reach the disk edge undamped. This incomplete damping is important for the 
question of the planetary migration. The migration speed and direction depend sensitively on a 
delicate balance between the amounts of the angular momentum which the planet deposits into the 
inner and outer parts of the disk. If the density waves dissipate completely and transfer all their 
angular momentum to the disk fluid then the only difference in torques acting on both sides of the 
disk is due to the surface density and temperature gradients in the disk, and to asymmetries in the 
locations of the inner and outer Lindblad resonances (Ward 1986). The amount of the resultant 
torque which leads to the orbital evolution of the planet is only ^ hp/vp compared to the magnitude 
of one-sided torque (GT80). Interaction with the outer part of the disk is usually stronger than 
with the inner one, leading to an inward migration (Ward 1986). 

Let us now assume that tidal perturbations are not damped completely upon reaching the disk 
edge and are reflected from it. The remaining waves will be dissipated on the way back to the 
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planet but some of them might survive and interact gravitationally with the planet, returning to it 
some of the initially launched angular momentum. If this effect is able to return to the planet about 
hp/vp of the one-sided angular momentum then the migration could be strongly modified [Tanaka 
&; Ward (2000) studied a similar effect caused by asymmetries in wave damping]. Consider, for 
instance, a planet sitting close to the outer edge of the disk, but still several hp from it (otherwise 
strong asymmetries in the torque will be produced when the wake is still forming). A tidal wave 
launched in the inner disk might be completely dissipated because of the large distance it has to 
travel in the disk. At the same time the wave in the outer disk might not shock at all before 
being reflected from the disk edge and can bring a significant amount of the angular momentum 
back to the planet. Thus, interaction with the outer part of the disk is now weaker than with 
the inner one and migration will change its direction - the planet will move outwards. In the 
same way one could show that planets near the inner edge of the disk tend to move inwards faster 
if the waves are incompletely damped in the inner part of the disk. One can roughly describe 
this process by saying that the planet tends to be pushed out from the disk towards its closest 
boundary. The same picture will hold for wave reflection off the edges of gaps formed by giant 
planets. These conclusions depend on a lot of assumptions, such as the details of the reflection 
process and gravitational interaction of the reflected wave with the planet, which certainly deserve 
further study. Whether this process is an interesting issue for the question of planet formation and 
survival in the course of migration depends on the relevant timescales. Nevertheless, incomplete 
damping of the density waves introduces additional degree of freedom on which the migration 
process depends and which could be important in some systems. 

Deposition of the wave angular momentum into the disk fluid leads to the evolution of the disk 
itself (Larson 1989). Spruit (1987) and Larson (1990) found that the action of shocked density waves 
is equivalent to an appreciable viscosity with corresponding dimensionless a-parameter (Shakura 
& Sunyaev 1973) reaching ~ 10~^ — 10~^. However they did not specify the source of the tidal 
perturbation, a deficiency remedied in GROl. There it was demonstrated that if all the solids in 
the disk were deposited into a population of Earth-sized objects then, again, an effective viscosity 
a ~ 10~^ — 10~^ is produced. 

Our results allow one to study the global evolution of the disk affected by all the planets 
present in it. The theory of the time-dependent accretion disks (Pringle 1981) states that the mass 
accretion rate at each point in the disk is uniquely determined by the divergence of the angular 
momentum flux carried by the waves: 



M 



Since fj depends on the distance travelled by the wave in a complex way [see equation (42)], one 
should expect M produced by a single planet to be a function of r. Using (42), (43), & (50) we can 
calculate the accretion rate in the disk at a distance r from the center produced by a single planet 
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located at r = r^: 



M{r) = sign(r - r,)j:^ (^J ^^(rp)So(r,)r^ (^^ j 
\ 3/2 V r \ P/"^-^/^ / M \ 



Here = d^{x)/dx. Note that because < 0, the accretion rate is positive inside of the 

planet (inflow) and negative outside of it (outflow) - the perturber tries to repel the surrounding 
gas. 

In Fig. 4 we plot the accretion rate in the disk due to the planetary torques from 8 planets of 
the Solar System (excluding Pluto) using equation (51). It is assumed that (Hayashi 1981) 

So(r) = 1700 g cm-^^rJl^ co(r) = 1.2 km s-^t'IJ^ (52) 

{rAU is the distance from the center measured in AU), and for the masses of the giant planets we 
take only the masses of their rocky cores: 15 for Jupiter and Saturn and 10 M© for Uranus 
and Neptune (otherwise they are likely to open gaps). 

One can see that in some locations in the nebula M ~ (1 — 5) x 10~^ Mq/jt suggesting 
a significant surface density evolution there. Indeed, the total disk mass inside Neptune's orbit 
calculated using equation (52) is only ~ 0.007 Mq. Averaging M over the bulk of the disk one 
obtains (M) ~ (2 — 4) x 10~^ Mq (depending on whether one averages M or \M\) which implies the 
typical dispersal time of the MMSN by planetary torques ~ (2 — 3) x 10^ yr, in rough agreement 
with observations (Hartmann et al. 1998). For more massive disks evolution could be more rapid, 
because the mass contained in the planets is increased (for a fixed disk metallicity) : planets could 
be more massive, meaning higher ratio Mp/Mi [see equation (51)], or simply be more numerous. 
This leads to stronger accretion so that the timescales < 10^ yr could be typical. Note however 
that at some point this increase in accretion rate is stopped by the tendency of massive planets 
to open gaps in the disk. This could drastically reduce their influence on the nebular evolution. 
Notice also that since M is very inhomogeneous radially, planetary-driven disk evolution must be 
highly time-dependent. 

A very useful feature of our analysis presented in §2.2 is that it can be applied to disks in 
which the surface density and sound speed have arbitrary radial distributions. They should only 
vary on scales larger than the perturbation wavelength for our tight-winding approximation to be 
valid (since in differentially rotating disks radial wavenumber rapidly grows with distance from 
the planet this condition does not pose serious restrictions). Thus one can not only calculate the 
instantaneous accretion rate at each point in the disk but also study the self-consistent temporal 
evolution of the disk driven by the planetary torques. This is very important, for example, for 
studying the gap formation around the planet in gaseous disks and we are going to investigate this 
process in a future study (Rafikov 2001). 
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Fig. 4. — Dependence of the planet-induced accretion rate M (in MQ/yr) upon tire distance r (in 
AU) in the minimum mass Solar Nebula (MMSN). Torques produced by 8 major planets (only 
masses of the rocky cores are taken for giants) are included here, and calculation is done using 
usual MMSN parameters [see equation (52)]. Positive values of M (inflow towards the center) are 
displayed in the upper panel, while negative ones (meaning the outflow) are in the lower panel. 
Dots denote the locations of the planets. One can see that in some positions in the nebula M could 
reach > 10~^ Mq/jt leading there to a significant surface density evolution on timescales ~ 10^ 
yr. 
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5. Conclusions. 

Tidal interaction of the gas disk with a planet embedded in it is important for the question 
of the orbital evolution of the planet as well as for the fate of the disk itself. In this paper we 
presented a global description of the evolution of density waves in vertically isothermal disks where 
two-dimensional fluid equations provide proper approximation. The disk surface density and sound 
speed are allowed to vary independently with radius. 

Our quantitative results are not very different from those obtained by Goodman & Rafikov 
(GROl): for low-mass planets surface density perturbations are weakly nonlinear and they shock 
after propagating several local disk scale lengths hp from the planet. Subsequent damping of the 
wave transfers its angular momentum to the disk and is intrinsically asymmetric, which could be 
important for planet migration. Disks evolve due to this angular momentum flux deposition and in 
the absence of other viscous mechanisms this could be the only driver of their evolution. We have 
demonstrated that for the parameters similar to those of the Solar System the tidal perturbations 
alone could produce spatially nonuniform and time dependent evolution with average accretion 
rates M ~ 10^^ — 10^*^ Mq/jt (yielding typical timcscale ~ 10® — lO'^ yr). In protoplanetary 
systems with more favorable conditions disk evolution may be even stronger. 

The prescription for the angular momentum deposition by planets in disks with arbitrary 
surface density and temperature profiles described here could easily be incorporated into other 
problems such as the gap formation around massive planet or planet-driven global evolution of the 
surface density in the disk. 
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A. Reduction of the equation for i?_ to Burgers equation. 

The equation for i?_ = u — 2c/ (7 — 1) = 2(co — 2c)/(7 — 1) reads 

[d^ + {u- c)d^] = - 



:Of:F — —O^Fq — co^ — cu-^ uo^- 



So ^ ^7-1 S ^7-1 

cu cv ^ ^ ^ c ^ , . 



We find it useful to introduce a new function 



7-1 Co 



Using this definition and equation (5) one can find that up to the second order in 

-V' + /"L V''- (A3) 



^-^0 = + 3-7 ^,2 

So 7 + 1 ^ ^ (7 + 1)2 



Also from the conservation of the Riemann invariant i?+ we have that 

u = 2^ = -2^^. (A4) 
7-1 7 + 1^ ^ ^ 

Now, we want to relate d^P to the derivatives of S — Sq- When doing this we have to remember 
that relation (4) holds true only in the reference frame comoving with the fluid. Thus, Eulerian 
increment of any quantity A^: has to be related to its Lagrangian counterpart Aj^hy Ae = A/,— dV, 
where d is a Lagrangian displacement. Obviously, V • d = — (S — So)/So) and since the disk is 
assumed to be thin and the spiral pattern is tightly wound, d^d !^ V-d = — (S — So)/So- Expanding 
Lagrangian increment of P up to the second order in S — Eq we obtain that 

d^P = d^Po + [c2(E - Eo)] + (S - So)^ + P^YT^^ §■ ^^^^ 

For a polytropic equation of state with a fixed polytropic constant (entropy in the disk is indepen- 
dent of r) the last term in (A5) is absent and the derivative of the pressure could be taken without 
worrying about the Lagrangian displacement. 

Since the tidal perturbation is assumed to be weak and tightly wrapped, the most impor- 
tant nonlinear terms responsible for the wake evolution are proportional to ipd^i/j while terms like 



This preprint was prepared with the AAS lAT^^X macros v5.0. 
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if)^ can be disregarded. Substituting (A3), (A4), and (A5) into (Al) we get after lengthy but 
straightforward calculation that 

- co(l + ip)d^ip 



+ 



7+1 
4co 

c 



■ n 

2 V H 

An AQr 



cu cv 



6c; 



ZCq 



7 + 



^V'S^lnco ^V^glnSo 

- 1 7 + 1 



(A6) 



We needed to expand S — Eq and d^P up to the second order in perturbation to make sure that 
all the terms proportional to d^'^ and ij^d^ip cancel out in the r.h.s of (Al). 



To an adequate approximation we can express from (26) 5<^v as 

2B 



An{r 



-u. 



(A7) 



Here we neglected the nonlinear terms ud^fyV and uv / {AVtr); the term with only slightly changes 
the propagation velocity of the wake and thus is disregarded too. 

In equation (26) the terms which are in phase with u or are important for the amplitude of 
the perturbation evolution while those out-of-phase only slightly affect the characteristic velocity 
but not the wave amplitude (as we noticed before). Derivative with respect to (j) changes phase 
by 7r/2; thus, terms with in equation (26) can be considered separately from others. We can 
integrate them over (p to obtain 
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1 



AQr 7+1 

This result, combined with equation (A4), confirms that v -^u (see §2.2). 
Using relations (A4), (A7), and (A8) we finally get that 



(A8) 



d4,ip - Co (1 + V') d^ip = Tp 



Co 
AQr 



1 1(9 In So sain Co A{r) 

2 2 91nr ^ 2 dlnr ~ An 



+ 0(^2). 



(A9) 



The local approximation of GROl could be retrieved now by assuming So = const and cq = const 
and expanding AO, to first order in terms of r — Tp. 

All effects of the nonlinear wake evolution arc embodied in the term tpd^tp in the l.h.s. of 
equation (A9). We now reduce this equation to the conventional Burger's equation following the 
approach outlined in GROl. First, we make a change of independent variables from ^ to 
given by the relations 



/ 



40 



(AlO) 



Here (j)' has the meaning of the azimuthal distance along the wake while r]i represents the displace- 
ment from the wake center in the ^-direction. Since the wake is narrow, rj <^ (j) and we can consider 
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co to be a function of ^' only: co(C) = co{4>). This transforms equation (A9) into (we drop the / 
from (f)') 

d^ip -ipdr,itp = -il;--^d^g{(f)), (All) 



where function g{(^) is defined by 



1 _^ l91nSo _^ 3 51nco ^(0) 



(A12) 



d4>X = X^dr,,x- (A13) 



2 2 51nr 2 91nr A0((^) 

Introducing new function x{4') = 9{4')'4^ we rewrite the equation (All) as 

1 

Changing from to a new variable t\ given by 
we arrive at the Burger's equation 

dtiX-xdrjiX = 0. 

Both ti and rji are dimensionless and we find it useful to rescale them to t, rj in the following 

way: 

V = ^fvi, t = %^, (A15) 
ip ip 

where Ip = Cp/\2A{rp)\ is a Mach-1 length at the position of the planet (in the limit r — >■ this 
rescaling makes our rj identical with rj used in GROl). Thus, we obtain equation (30). 

Finally, we calculate the behavior of function g. From definitions (24) and (AlO) one finds that 
d(f) = — (Ar2/co)(ir. Using this and the definition of Oort's parameter A one can easily integrate 
equation (A12) to find 

5 = Ci(^^) ' ^"1=0071.^. (A16) 

The arbitrary constant Ci should be chosen in such a way that in the limit r ^ Vp our definitions 
of X and t reduce to the analogous expressions of GROl [given by equation (35)]. One achieves this 
by taking 

21/4 

Ci = (A17) 



[see equations (34)-(35)]. 



